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ABSTRACT 

We consider the problem of convergence in homogeneous shearing box simula- 
tions of magneto-rotationally driven turbulence. When there is no mean magnetic 
flux, if the equations are non dimensionalized with respect to the diffusive scale, 
the only free parameter in the problem is the size of the computational domain. 
The problem of convergence then relates to the asymptotic form of the solutions 
as the computational box size becomes large. By using a numerical code with a 
high order of accuracy we show that the solutions become asymptotically inde- 
pendent of domain size. We also show that cases with weak magnetic flux join 
smoothly to the zero flux cases as the flux vanishes. These results are consistent 
with the operation of a subcritical small-scale dynamo driving the turbulence. 
We conclude that for this type of turbulence the angular momentum transport 
is a proportional to the diffusive flux and therefore has limited relevance in as- 
trophysical situations. 

Subject headings: accretion disc - MRI - MHD - dynamos - turbulence 

1. Introduction 

The magneto-rotational instability (MRI ) is commonly invoked to explain the origin 



of turbulence in electrically conducting discs fjBalbus & Hawleylll99l[ ). Although the insta- 



bility itself is amenable to an analytic treatment, much of what is currently known about 
the nonlinear development of this instability is based on numerical simulations. In partic- 
ular, the efficiency of MRI driven turbulence at transporting angular momentum, which is 
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the crucial quantity that ultimately controls the accretion rate, is known almost exclusively 
through numerical simulations. Because of the difficulties inherent in simulating an entire 
disc the majority of the num erical work has bee n based on the local model known as the 
shearing-box approximation (IHawley et al.l 119951 ) . The advantages of the shearing-box ap- 
proximation (SBA) are numerous, the cylindrical geometry of the full disc is replaced with 
the Cartesian geometry of a rectangle, simple periodic and shear-periodic boundary condi- 
tions can be applied, and most importantly, whatever numerical resolution is available can 
be deployed to resolving the ensuing turbulence. However, a number of issues arise that have 
led several authors to question the applicability o f the s hearing-box approximation (SBA) 
to the study disc turbulence. iRegev &: Umurhanl (120081 ) have noted some inconsistencies 
in the formulations of the SBA with uniform background states, they have also questioned 
the assumption of locality t hat underlies the d erivation of the SBA. Related concerns about 
locality were expressed by iBodo et al.l (120081 ). However, the biggest issue about the SBA 
remains the so-called problem of convergence. 

There are two configurations in which MRI driven turbulence is studied numerically; one 
in which there is a net magnetic flux threading the layer, the other in which there is none. 
In the former case, if the uniform component of the magnetic field is vertical, say, there 
is a linearly unstable mode with a well defined vertical wavenumber of maximum growth 
rate that sets the scale of the instability, in the latter no such state exists, and the MRI 
must set in as a subcritical instability. Similar considerations apply if the initial field is 
azimuthal. What exactly determines the characteristic scale of the turbulence in the no- 
net-flux case is, at the moment, an open question. The conceptual appeal of the no-flux 
case is that it offers the possibility of a universal state of MRI turbulence in which the disc 
becomes self magnetized through dynamo action, so that the angular momentum transport 
depends on the disc properties but not on the amount of flux threading the disc. Clearly 
much effort has been devoted to determining if such a universal state exists, and the value 
of the associated turbulent transport. Simply stated, the problem of convergence is that the 
angular momentum transport measured in numerical simulations based on the SBA with 
homogeneous background state and zero mean-flux depends on numerical resolution, and 
decreases as the resolution increases. 



This effect was first noted by iFromang fc Papaloizoul (120071 ) in a series of numerical 
experime nts with the ZEUS code and subsequent l y confirmed by ot her authors using different 
codes (eg iPessah et al.l 120071 ; ISimon et al.ll2009l ; iGuan et al.ll2009l ). There are a number of 
important points that should be made. Most of the evidence for the convergence problem 
is based on codes with no explicit viscosity or magnetic diffusivity in which the dissipation 
arises solely from truncation errors. An increase in resolution is interpreted as a decrease in 
effective dissipation, so that the convergence problem can equivalently be stated as a decrease 
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in the effective transport with decreasing dissipation. Recently, however, iFromangj (120 101 ) 
indicated that the problem does not arise when explicit viscosity and magnetic diffusivity are 
included. Also, the problem is seen in simulations with ho mogeneous background s tate and 



periodic boundary conditions in the vertical. Simulations by lKapyla fc Korpil (120101 ) suggest 
that the convergence problem is absent if different boundary conditions are applied in which 
the magnetic field is vertical on the upper and lower bounda ries. Also, the problem seems to 



disappear when stratification is included (IDavis et al.ll2010l ). or the aspect ratio is such that 



that computational domain consists of a sufficiently tall box (Stone, private communication). 
Finally, even when the effect is clearly observed there is no universal agreement about the 
rate at which the transport decreases with resolution. A number of questions naturally arise: 
what is the origin of the convergence problem? Is it a numerical artifact associated with the 
indiscriminate use of ideal solvers? Is it a real physical phenomenon related to subcritial 
dynamo instabilities? Is it that the shearing-box approximation in its simplest form is too 
idealized to describe MRI turbulence in a disc? What is the asymptotic scaling of the angular 
momentum transport with resolution/dissipation, and what can we learn from it? 

In this paper we address some of these issues. We revisit numerically the problem of MRI 
driven turbulence in shearing-boxes with uniform background state and periodic boundary 
conditions in the vertical. We formulate the problem in a slightly unconventional way that 
is designed to emphasize the role of symmetries in the shearing-box approximation, and to 
make it easier to distinguish between numerical and physical effects. 

In the ideal limit, the shearing-box equations are scale invariant, in the sense that dou- 
bling the number of grid points in a simulation, or doubling the size of the computational 
domain are entirely equivalent. This gives rise to two equivalent limiting procedures; letting 
the grid spacing vanish for finite domain size, and letting the domain size become infinite 
for finite grid-spacing. Although the two limits are equivalent, the latter makes it easier 
to interpret the results in terms of some classical results from dynamo theory. Using this 
formulation we consider the results of high-resolution numerical experiments and reach the 
following conclusions. The convergence problem is real in the sense that the transport by 
MRI driven turbulence does decrease with decreasing dissipation and vanishes for vanishing 
dissipation. The problem is related to two aspects of the shearing-box approximation with 
periodic boundary conditions, one is the inherent symmetry of the shearing-box approxi- 
mation , the other is the absence of an effective inverse cascade in the dynamo solutions. 
We find that there exists an asymptotic regime in which the effects of the symmetries are 
apparent but it is realized at very high resolution, or as we shall see presently, at very large 
system size, suggesting that most of the simulations described in the current literature are 
not in that regime. Finally, we argue that even though the shearing-box approximation in its 
simplest form is probably unable to capture the physics of MRI driven turbulence in discs, 
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it nevertheless gives us very useful clues of which effects are likely to be important, and how 
to proceed to improve our models. 

The paper is organized as follows: in the next section we discuss the zero-flux case; we 
review the standard formulations and propose a slightly different approach that better sepa- 
rates between physical and numerical quantities. We present numerical evidence that suggest 
the existence of an asymptotic state in which the transport becomes independent of system 
size. We then discuss cases with small but finite flux, and show that they match smoothly 
with the zero-flux cases. Finally, we discuss the implications for modeling astrophysical discs. 



2. The zero net flux case 

We begin by discussing the case with no mean flux. Conceptually, this is the simplest 
case since for an ideal incompressible fluid, the SB A is completely scale invariant. As we 
shall see presently, this symmetry plays an important role in determining the asymptotic 
form of the solutions. 



2.1. Formulation 



The first step is to cast the SBA equations in dimensionless form. Often, this is accom- 
plished by selecting units that depend on the sound speed, although this choice is convenient 
from an astrophysical point of view it is not the most suitable to bring out the natural 
symmetries of the equations. Instead, we begin by considering an incompressible fluid-i.e. 
a fluid with infinite sound speed, with finite viscosity and magnetic diffusivity, and we shall 
re turn to the compress ible case p resently. A detailed prese ntation of the SBA can be found 



in 



Hawley et al.l (Il995[ );(see also iLesur & Longarettil 120071 ). The equations in dimensional 



form can be written as 



dv 

~dt 



v ■ Vv + 2Q x v 



Anp p 



B- 



•>7T 



+ P ) - V (2Attx 2 ) + z/V 2 v, 



(1) 
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V-v = V-B = (3) 
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where B, v, p and P denote, respectively, magnetic field, velocity, density and pressure. The 
local angular velocity fi = Qe z and shear rate 

A = -¥ (4) 
2 OR K J 

are considered constants. For a Keplerian disk f2 oc R~ 3 / 2 and A = — (3/ 4)0. The velocity 
can be decomposed as the sum of the base Keplerian flow and the fluctuations 

v = 2Axe y + u, (5) 

likewise, the pressure can be decomposed as the sum of the average and the fluctuations 

P = Po+P- (6) 

With these decompositions, the dimensional shearing box equations become 

Du * 9u . n „ B-VB 1^/B 2 \ - 

— + Ax— + Au x e y + 2Qxu = — V — + p + z/V 2 u, 7 

Dt ay Anp p \8n J 

Z)B <9B 

+ Ax— - AB x e y - B Vu = V V 2 B, (8) 
Dt ay 

Vu = VB = 0, (9) 

with D/ Dt = d/dt + u • V. The above system of equations contains only three dimensional 
parameters: the rotation frequency Q, the viscosity v and the magnetic diffusivity r\. 

By adopting the following units of time, length, velocity and magnetic field intensity 
r = ^; C = l D = u* = v 7 ^; B* = y/puQ. (10) 

the equations can be written in dimensionless form as 

Du 3„<9u 3. n A 1 A - a/b 2 A ~„ 

— r - -x— - -u^ + 2e,xu = — B • VB - V — + p + V 2 u 11 

Dt 4: dy 4 y 4tt I 8tt / 

L>B 3„<9B 3 A - „ 1 , . 

7rf-4^ + 4 B ^- B ' Vu+ ^ VB = ° < 12 » 
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where the hatted quantities are dimensionless and P m = v/rjis the magnetic Prandtl number. 
It is important to note that the induction equation (j!2p only depends on P m , while the 
momentum equation ffTTl) has no adjustable parameters at all. This is because in the SB A 
with no mean flux, the diffusive scale Id is the only intrinsic length that can be defined using 
the dimensional parameters of the problem. Clearly, in order to solve equations ffTTl - TT21) 
one also needs to define a computational domain and suitable boundary conditions thereby 
introducing an external length-scale L that characterizes the size of the computational box; 
the latter, however, has no direct physical meaning. In particular, in the case of periodic 
boundary conditions (shearing periodic in the radial direction), formally one is considering 
an infinite domain, with L representing a long wavelength cutoff. The domain size can itself 
be nondimensionalized to define the additional parameter R = L/C. For fixed magnetic 
Prandtl number, R- the domain size in units of the diffusive scale -is the only adjustable 
parameter. 

We note that the non-dimensional form of the equations above is different from what 
is typically found in the literature. In general, the box size is used as unit of length and 
the square of the parameter R is then identified with the Reynolds number. However, in a 
turbulent flow, the Reynolds number should be more properly defined in terms of the integral 
scale of the turbulence. The question then becomes whether asymptotically, for R — > oo, 
in MRI turbulence the integral scale is linked to the box size, or to the diffusive scale, or 
to a combination of both. We note that if the solutions become asymptotically independent 
of the box size, the system of equations ffTTl - fT21) do not depend on any non-dimensional 
parameters, apart from the Prandtl number, i.e. for fixed Prandtl number there is a universal 
solution. 

In the zero net flux case, dynamo processes must effectively maintain the magnetic field 
against dissipation, however the magnetic field is itself responsible for generating, through 
the MRI, the turbulence that sustains it. It is well known that, in a (small-scale) turbulent 
dynamo, the magnetic energy is concentrated near the resistive scale. In this case, one 
expects that the MRI driven velocity may also be characterized by the diffusive scale. In 
order for the velocity to have a significant component comparable to the box size an efficient 
mechanism is required that can transfer energy backwards, from the resistive scale to the 
largest available scale, i.e. an inverse cascade. Whether such an inverse cascade exists is at 
the heart of the problem. 

As an aside, we note that the condition for the validity of the SBA is that the char- 
acteristic scales of the solution be much smaller than the box size, lest the use of periodic 
boundary conditions is not justified. Therefore, if there were cases in which the shearing box 
results would show a dependence of the solutions on the box size R as R — > oo, one should 
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actually reformulate these cases by introducing some global scale characteristic of the whole 
disc. 

One of the main objectives of MRI studies is to determine the efficiency of the transport 
of angular momentum by Maxwell and Reynolds stresses. It is customary to define the 
quantity 

„ B T Bi, 

E = < u x u v > 13 

that represents the box and time averaged value of the total stresses (hereinafter overbars 
denote time averages, while angle brackets denote a box average). This quantity has di- 
mensions of square velocity, thus with our choice of non-dimensionalisation it is measured in 
units of VLv. For fixed magnetic Prandtl number we can write 

s ~ f(R)n% = f(R)nu. (14) 

Here, f(R) is a dimensionless function that accounts for the possible dependence of the 
solutions on the externally imposed box size. When R = 0(1) the box size is comparable 
to the diffusive scale, and, on general grounds, one expect a strong dependence of f(R) on 
R. It is not clear, a priori, what dependence one should expect as R — > oo. It is important, 
however, to be clear about the relationship between the asymptotic form of f(R) and the 
transport efficiency. In turbulence theory, it is customary to refer to "turbulent transport" 
to describe a transport process that becomes independent of diffusion for small diffusion, 
and to "collisional transport" to describe one such process that remains proportional to 
the diffusivities when these become small. Because, here, £ is measured in diffusive units 
(fiz/), if f(R) ~ const, as R — > oo the angular momentum transport is entirely collisional. 
Conversely, in order for the transport to be "turbulent" in the sense defined above, f(R) ~ R 2 
as R — > oo. In fact, any asymptotic dependence of f(R) weaker than R 2 , will ultimately 
lead to a transport that vanishes for vanishing viscosity. 

Typically, in the literature, the transport efficiency is measured in terms of the parameter 
a. In the incompressible c ase this is defined as the total stresses measured in units of L 2 Q 2 



(ILesur fc Longarettil 120071 ); in other words 



— ^ = f( R >w (15) 

since R 2 = L 2 Q/u. Again, we recover the result that in order for a to have asymptotically 
finite (constant) value the function f(R) must diverge quadratically as R becomes large. 

When compressibility is taken into account, an additional physical parameter is intro- 
duced, namely the sound speed c s , that can be used to define a new length scale H = Qc s . 
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When vertical gravity is accounted for, H, represents the scale height of the disc. In stud- 
ies such as the present one, in which vertical gravity is neglected, this length has no direct 
meaning, nevertheless it introduces a new non-dimensional parameter Rh = H/C, that must 
be taken into account in the scaling arguments discussed above. As before we can write 

E = F(R, R H )Qu, (16) 

where F(R,Rh) is a dimensionless function. For Rh = 0(1), again we expect a strong 
dependence of F on Rh- However for Rh large, F must become asymptotically independent 
of Rh, since this corresponds to the incompressible limit in which F must approach /. 

Physically, this corresponds to situations in which H is much larger than the character- 
istic scales of the solutions. Again, if there is no dependence of these scales on the externally 
imposed box size L, the only available length is £ and the characteristic scales of the solu- 
tions will be proportional to C. For constant sound speed and therefore constant ratio H/L, 
the limit R — > oo corresponds to Rh — > oo, and therefore the st ress behavior should be the 



same as in the incompressible case. In this case a is defined as in lShakura Sunyaevl (119731 ) 



and we therefore expect again the scaling 

1 



OL ~ . 

R 2 



2.2. Numerical Considerations 

In most of the MRI studies, as well as in this paper, the analysis of MRI turbulence is 
performed using ideal codes, where dissipative effects only arise due to numerical truncation. 
In general, in this case, it is not known how to write the dissipation term explicitly, however 
it can be said that the amount of dissipation depends on the ratio between the local scale of 
variation of the physical quantities and the cell size. Structures comparable to the cell size 
are subject to strong dissipation that decreases when the scale of the structures increases. 
Different numerical schemes, however, will behave in different ways, in particular the reduc- 
tion of dissipation when the local scale is decreased will be sharper for high order schemes 
and gentler for lower order ones. In (fTUj) we defined the dissipation scale in terms of the 
viscosity and the shear as Id = C = a/ u/Q. In the numerical approach one is not able to get 
an explicit expression for the dissipation scale, nevertheless it is known that Id is related to 
the cell size 5, in a way that, however, depends on the scheme. In addition, it is know that 
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Id (in units of 5) will be larger for lower order schemes and smaller for the high order ones. 
Thus, in the non-dimensionalization of the equations, in principle one should use Id, but in 
practice one can only use 5. With the two dimensional parameters Q and S, it is possible to 
define the units of time r, length C, velocity u* and magnetic field B* as 

r = -; C = 5; u* = Q5; B* = ^/pQd (19) 

By the same considerations discussed above, the scaling of stresses can be written as 

S ~ f{N)Q 2 5 2 (20) 

where N = L/5 represents the number of grid points and the f(N) factor accounts for a 
possible dependence of the solutions on the externally imposed box size. The scaling for a 
is correspondingly given by 

a ~f( N )]p- ( 21 ) 

If there is no dependence on the box size, i.e. f(N) is constant, we have aN 2 = const; this 
constant however will be different from scheme to scheme, since different numerical schemes, 
as discussed above, will have different relationships between the dissipation length l D and 
the cell size S. The proper unit for S would be Q 2 1 2 d, therefore we expect 



aN 2 ~ (J) . (22) 

Since Id is expected to be smaller for high order codes the same should be true for the value 
of aN 2 . It should be noted that this discussion ignores possible difference in the effective 
Prandtl number for the different schemes, that can account for additional residual differences 
between codes. 

For the compressible case, we can repeat the same reasoning as above: we have an 
additional parameter Nh = H/S that represents the number of grid points over H , however 
the dependence on this parameter should disappear, for constant sound speed, in the limit 
iV ->■ oo. 



7 x 2 
Li 



2.3. Numerical results 



We have performed a series of compressible isothermal shearing box simulations for 
the zero net flux case, with di fferent resolutions an d with three different numerical schemes 
available in the PLUTO code (IMignone et al.l 120071 ) . The main difference between the three 
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schemes is the order of accuracy of the reconstruction, more precisely we employed a scheme 
with TVD linear reconstruction (second order accurate for smooth flows), a parabolic re- 
construction (PPM) and a monot onicity preserving fifth o rder reconstruction (MP5) (fo r 
a description o f this scheme see ISuresh &: Huynhl (119971 ); iMignone fc Tzeferacod ( 120101 ); 
Mignone et al.l (120101 ) ). Our goal is to study the asymptotic behavior of the solutions as 
the separation of scales between the box size and the dissipation length becomes large. This 
can be achieved by increasing the resolution for a fixed scheme, or for a fixed resolution, by 
using a scheme with higher order of accuracy. 

In the shearing box approach, the Cartesian coordinates x, y and z refer respectively 
to the radial, azimuthal and vert ical directions. Our computat ional box has aspect ratio 



L x : L y : L z — 1 : 7r : 1. As in iFromang &: Papaloizoul ( 120071 ) we use cells elongated in 



the azimuthal direction, with aspect ratio 1:2:1. For each of the schemes we used four 
different resolutions: 32, 64, 128 and 256 points in the vertical direction, the corresponding 
grid sizes are therefore 32 x 50 x 32, 64 x 100 x 64, 128 x 200 x 128 and 256 x 400 x 256. 
The sound speed is chosen such that H = L, the initial magnetic field is (0, 0, B sin(27nr)), 
with B corresponding to = 1500 and random noise in the y component of the velocity is 
introduced initially to start the growth of the instability. 

As it is well known, the time histories of quantities like the box averaged magnetic 
energy or the box averaged Maxwell or Reynolds stresses, have an initial transient followed 
by the establishment of a stationary state in which the quantities fluctuate around a well 
defined mean value. The amplitude of the fluctuations decreases with increasing resolution. 
To estimate the mean value we average over the simulation time excluding the initial tran- 
sient phase. The simulation time varies from case to case and it is always larger than 100 
revolutions. We estimate the error of the mean value by subdividing the averaging interval 
in ten subintervals and then computing the standard deviation of the subinterval averages. 
The resulting estimate of the error is always less than 10%. 



2.3.1. Transport 

We start our analysis by considering the behavior of a, as defined in (1171) . as a function 
of resolution. Figure [1] shows aN 2 as a function of N for all the cases considered. The 
quantity displayed represents the function f(N) defined in (}2"Tj) . which accounts for a pos- 
sible dependence of the solutions on the externally imposed box size. Consonant with our 
objectives, we are interested in the asymptotic behavior of f(N) as iV — > oo. 

We recall that in order for a to be resolution independent, f(N) should diverge as N 2 . 
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The results by iFromang fc Papaloizoul (120071 ) suggested that f(N) oc N. A constant value 



of f(N), on the other hand, means that the solution is asymptotically independent of the 
box size and that the scaling of a is determined solely by the equations. 

In Figured! the squares refer to the linear scheme, the stars to the PPM and the triangles 
to the MP5. It is obvious that, as the the order of accuracy of the scheme is increased, the 
curves tend to be flatter for increasing N. However, while the linear and PPM schemes show 
a continually decreasing slope and do not seem to have reached yet an asymptotic behavior, 
the MP5 scheme, by contrast, appears to be leveling out. The results of MP5 suggest that 
f(N) approaches a constant value as iV — > oo and therefore that a has the asymptotic 
natural scaling ~ 1/N 2 arising from the equations alone. Our expectation is that, given 
sufficient resolution, a similar behaviour would be observed with the other schemes. 

It is important to note that even if all the curves corresponding to different schemes 
eventually level off, they will not in general have identical asymptotic values. The reason for 
this is that we do not know the explicit form of the numerical dissipation for each scheme; 
accordingly we have used the cell size instead of the dissipation length to non-dimensionalize 
the equations. Indeed, we note that in agreement with (1221) the higher order scheme MP5 
has a lower value of the non-dimensional stresses than the other schemes. 



2. 3. 2. Structures 

The behavior of a as a function of resolution suggests that the only scale that determines 
the properties of the solutions is the dissipative scale, which, in numerical studies, is related 
to the cell size. It is therefore important to examine the behavior of the size of the typical 
observed structures as a function of resolution. In order to characterize the scales of magnetic 
structures, we measure the average scales of variation of the field in the directions parallel 
and perpendicular to itself. We define the quantities 



1/2 / < | B |4> y/2 
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which represent measures of the characteristic scales of the magnetic field respectively in 



the p arallel direction and along the direction of maximum gradient (see ISchekochihin et al. 



20041 ) . In fig. [2] and [3] we plot respectively Nl± and iV7|| as functions of N. As before, 
the three curves refer to the three different numerical schemes and the symbols have the 
same meaning as in the previous figure. The lengths are normalized with respect to the 
cell size and the plotted quantities therefore represent the average number of grid points 
on which the observed structures extend. As expected, the magnetic structures are highly 
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Fig. 1. — Plot of aiV 2 as a function of iV. The three curves refer to the three different 
numerical schemes used in the simulations. These are identified by the symbols in the 
legend. In all cases, the error bars are smaller than the sym bols. The slope is about 1, 
corresponding to the results by iFromang fc Papaloizoul (120071 ) with the ZEUS code, for the 
PPM curve between N=64 and N=128 . 
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anisotropic with l\\/l± ~ 10. Considering the dependence on resolution, we see again that the 
linear and PPM schemes show an increase of the number of grid points over which magnetic 
structures extend, while the MP5 scheme appears to tend to a constant number of grid points 
independent of resolution, both in the longitudinal and transverse directions. 



100 



10 



O linear 


i 




MZ PPM 






- A MP5 




























1----A 

i 


A - 



100 
N 



Fig. 2. — Plot of A/7 1 1 as a function of A/. The three curves refer to the three different 
numerical schemes used in the simulations. These are identified by the symbols in the 
legend. 

This analysis suggests that the magnetic field is characterized by elongated structures 
extending transversally over a few grid points. In order to get further insight into the three- 
dimensional structure one can define a third scale in the direction orthogonal to both B and 
B x J and that points approximately in the vertical direction: 



1/2 

*B-J=( TTp Tl2 \ ) • ( 24 ) 
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Fig. 3. — Plot of iV7_|_ as a function of N. The three curves refer to the three different 
numerical schemes used in the simulations. These are identified by the symbols in the 
legend. 
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The ratios of the scales along the three directions are tipically l± : lu : Zb j ~ 1 : 10 : 5 
indicating th at the magn e tic st ructures can be thought of as thin sheets. This result is also 



discussed by iGuan et al.l ( 120091 ). who base their analysis on the correlation function. They 
find that the correlation lengths scale as iV -2 / 3 and that therefore, measured in units of 
the cell size, show an increase with resolution, suggesting that one could expect some kind 
of transition. We have also repeated this convergence analysis on the correlation lengths 
and again we find that while the linear and PPM schemes show a behavior similar to that 
discussed by Guan et al. (2009), the MP5 scheme at a resolution of 256 grid points shows a 
convergent behavior, i.e. the correlation lenghts scale as the cell size. A better impression 
of the magnetic field structure can be obtained by looking at Fig. HI where we show a 2D 
cut of the Maxwell stress distribution in the x — y plane at z — . The 3D sheets appear in 
the figure as filaments of high Maxwell stresses and high magnetic field intensity, separated 
by regions of low magnetic field. 

We can now discuss this filamentary structure from a more physical point of view. 
From the y-component of the induction equation, it can be seen that the field structure is 
determined by the competition between three terms: the stretching by the background shear, 
the stretching due to velocity fluctuations (B • Vv term) and dissipation, and the transverse 
size is determined by the balance between these terms. The background shear tends to 
stretch the field along the azimuthal direction, simultaneously decreasing its perpendicular 
scale of variation, until dissipation sets in. These magnetic filaments, however, give rise to a 
local transport of momentum and, therefore, tend to reduce the local shear. One, therefore 
expects that the effect of the nonlinear term B • Vv is to reduce the stretching by the 
background shear. The observed increase of the normalized transverse scale with increasing 
resolution may be related to the increase of this effect, which however, ultimately, has to be 
limited since it cannot completely suppress the background shear. Thus, asymptotically the 
transverse scale has to remain constant and be determined by the balance between shear 
stretching and dissipation. The MP5 scheme, which is closer to the asymptotic behavior, 
accordingly shows a tendency of the transverse scale to become constant. This leads to the 
conclusion that a possible dependence on the scale of the box may appear only in the parallel 
scale. However, the ratio l\\/l± seems to be independent from resolution, i.e. the magnetic 
structures do not appear to become more elongated as the resolution increases and all the 
characteristic lengths appear to become asymptotically constant when measured in unit of 
the dissipative length. 
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Fig. 4. — 2D cut in the x — y plane of the distribution of Maxwell stresses. The distribution 
is characterized by thin filaments with high stress values separated by wide regions where 
the stresses are almost zero. 
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2.3.3. Probability Density Functions 

Further details on the solutions can be examined by considering the behavior of the 
probability distribution function (pdf) of magnetic field intensity and of the second order 
structure functions. Our purpose is twofold: on one hand we want to show that, indeed, 
if we measure all quantities in their proper units, there is a convergence to a well defined 
solution when we increase the number of grid points and, on the other hand, we want to 
better characterize this solution. We start by examining the probability distribution function 
(pdf) of magnetic field intensity. Figure [5] compares the two MP5 cases with the highest 
resolution, the pdf's are normalized to unity and the field values are measured in units of 
B*. 

The pdf's have an exponential character, the mean value of |B|/S* is about 5.5 with a 
difference of less than 1% between the two cases at different resolution. The contribution to 
the magnetic energy and to the Maxwell stresses comes mainly from the high field values; 
values to the right of the vertical line plotted in the figure contribute to about 70% of the 
total magnetic energy and to about 80% of total Maxwell stresses, while occupying only 
about 20% of the volume. The higher contribution to the Maxwell stresses with respect to 
magnetic energy is explained by the fact that these higher field values correspond to the 
magnetic structures discussed above. For these, there is a strong correlation between the 
radial and azimuthal components of the field, while for lower field values the correlation 
is much lower. The fact that the volume fraction supporting the transport is independent 
of resolution explains the decrease in the amplitude of the fluctuations when resolution is 
increased. In fact, as resolution is increased and the volume fraction remains the same, the 
number of points that contributes to the average gets larger and therefore the fluctuations 
of the average become smaller. 

In fig [6] we plot the time averaged second order longitudinal and transverse structure 
functions for the y component of magnetic field, for the two MP5 cases at the highest 
resolution. They are explicitly defined as 



S 2 i = < (B y (x, y + h,z)~ B y (x, y, z)f > S 2t = < (B y (x + h, y, z) - B y (x, y, z)) 2 > 

(25) 

Following our discussion, we compute S 2 i and S 2t in units of B* 2 and h in units of 5. Again, it 
can be seen that, using these units, the difference between the two cases at different resolution 
is very small, indicating a convergent behavior. The structure functions in both directions 
become flat at scales larger than few grid points. This plateau indicates that the values of 
magnetic fields become uncorrelated at these scales. A simil ar conclusion could be reache d 



from the behavior of the averaged spectra, like those shown in lFromang &: Papaloizoul ( 120071 ). 
that become flat at small wavenumbers. 
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Fig. 5. — Probability distribution functions of the modulus of the magnetic field measured in 
units of B* for two MP5 cases with different resolutions. The number of grid-points per unit 
length-measured in units of L-are 128 (solid line) and 256 (dashed line). The field values 
larger than that indicated by the vertical line contribute to about 70% of magnetic energy 
and 80% of Maxwell stresses. 
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From these two results, we conclude that the solution is characterized by the super- 
position of uncorrelated magnetic sheets. Those with high field intensity also have a high 
degree of correlation between magnetic field components, and consequently give the highest 
contribution to both the total magnetic energy and Maxwell stresses. An increase of the do- 
main size does not change the properties of these sheets, but simply increases the available 
volume and therefore the number of sheets, leaving unchanged all the average properties. In 
a similar way there should not be any modification of the solution by changing the aspect 
ratio of the computational box, and we have indeed tested this by examining two further 
cases with aspect ratios respectively 1 : 2ir : 1 and 2 : ir : 1, with MP5 and 128 points in the 
vertical direction, finding no changes in all the averaged properties. 



3. The finite net flux case 



We now turn to the case in which a net magnetic flux threads the computational domain. 
In this case the magnetic field has a uniform component of strength B Q , say, that can be 
used to define a new dimensional length scale Xa given by 



A 



1 Bn 



(26) 



Accordingly, the averaged stresses will show a dependence also on this scale and keeping the 
same units as in the previous section, in general we can write 



9 



Id' 



D 



(27) 



Here, we are mostly interested in the form of the function g, when L/ljj — R — > oo. On 
general grounds, we expect, it will take different forms for different relationships between 
the three scales L, Xa and Id- We start our discussion by considering what is known from 
numerical simulations in the available literature. The case that has received the greatest 
attent ion is that of a mean vertical field; the results have been summarized in iPessah et al. 
(120071 ). For this case the system is linearly unstable to the MRI; the (vertical) wavelength 
of maximum growth rate is given by A m = 27Ta/16/15A^, the instability exists only for 
wavelengths A > 1/\/3Xa- If Xa is larger than a/3L, no u nstable wavelengths fit in the box, 
the system becomes stable, and the stresses drop to zero. IPessah et al.l ( 120071 ) show that as 
Xa increases, the stresses scale linearly with A^, reach a maximum and, eventually, drop to 
zero when the above condition is met. Furthermore, it has been f ound that, at fixed Xa/L 



that the stresses increase with resolution — i.e. with R (see e.g. ISilversI 120081 ; iBodo et al. 



20081 ); thus, for high enough resolution, they should converge to a value independent of R. 
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Fig. 6. — Longitudinal (left panel) and transverse (right panel), time averaged, second 
order structure functions of B y . The two curves correspond to two MP5 cases with different 
resolutions. The number of grid-points per unit length-measured in units of L-are 128 (solid 
line) and 256 (dashed line) 
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Also, we have to note that iBodo et al.l (120081 ) have found a dependence of the solutions on 
the aspect ratio L x /L z of the computational box, with a convergent behavior for high enough 
values of L x /L z . Putting these considerations together, we conclude that for l D <C A^ < L, 



^ ) n 2 L 2 



(28) 



where the function g is initially proportional to Xa/L, then reaches a maximum (at a fixed 
value of Xa/L) and finally drops to zero for A^ > \[ZL. Thus the value of E at the maximum 
scales as Q 2 L 2 . 

In the opposite regime, when I n ~ Xa « L, the averaged stresses tend towards the 
value obtained in the zero flux case (jPessah et al.ll2007l ). The results of the previous section 
suggest that, for large R, S is independent of L and scales as VL 2 l 2 D . 

In order to connect smoothly the two portions with, respectively, In ~ A^ <C L and 
Id « Xa < L, there must be an intermediate range of values of A^ in which the stresses 
scale as Q 2 X 2 A . In this range the solution becomes insensitive both to Id and to L, this 
can be achieved only if Aa is sufficiently far from both scales and therefore the extent of 
this intermediate range with quadratic scaling should increase as R increases. Contrariwise, 
decreasing R should cause the interval to shrink eventually to disappear for low values of R. 
In fig. [7J we show the results of calculations with a net vertical flux performed at a resolution 
of 128 points in the vertical direction, with the MP5 scheme. The horizontal dashed line 
corresponds to the averaged stress value for the zero flux case, the dashed line corresponds 
to a quadratic scaling, the solid line corresponds to a linear scaling and the vertical solid 
line marks the stability boundary; to the right of it, the stresses vanish. 

It is apparent from the figure, that indeed, for small enough Xa, a range of values 
of exists with quadratic scaling, while for larger A^ the curve reverts the to linear scaling 
discussed above. 

We have deliberately restricted our investigation to values of A^ for which we have 
quadratic scaling and the transition from quadratic to linear. 

Larger values of A^, for which the stresses increase linearly up to a maximum and drop 
to zero to the right o f the stability boun dary have already been extensively discussed in 
the literature (see e.g. iPessah et al.l 120071 ). We note here, that the portion with qu a dratic 
scaling has never before been reported in the literature (however lLongaretti fc Lesurl (120101 ) 
find evidence of a scaling steeper than linear), this is most likely due to limited resolution 
since, as discussed above, for low values of R, the quadratic portion disappears. We believe 
that the reason why we were able to detect it is because of the high order of accuracy of the 
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MP5 scheme, which for the same resolution, produces a much higher effective value of R. 



4. Discussion and Conclusion 

In this paper, we have addressed the problem of convergence in numerical studies of 
MRI driven turbulence in the homogeneous shearing box approximation. Any simulation 
that is physically useful should have some meaningful property that becomes independent 
of the number of grid points as the latter becomes large. This does not seem to be the 
case here. As noted by several authors, in the weak flux regime-i.e. when the total mag- 
netic flux threading the computational box is small or zero, the dimensional angular mo- 
mentu m transport decreases as as the resolution increases eventually becoming arbitraril y 



small (IFromang fc Papaloizoull2007t iPessah et al.ll2007l ; ISimon et al.ll2009t iGuan et al.ll2009l ) . 



There are at least two possible frameworks in which to address this issue: one in terms of 
Reynolds number, the other in terms of computational domain size. In both cases one is 
interested in the asymptotic behavior of the solutions as either the Reynolds number or the 
domain size become large. Clearly, both frameworks are related and to some extent it is 
a matter of preference which one it is picked. In this paper we have chosen to discuss the 
convergence problem in terms of domain size. The reason is that, because of the symmetries 
of the homogeneous shearing box approximation, there is no a priori characteristic outer 
scale for the velocity, which leads to an ambiguity in the definition of the Reynolds number. 
On the other hand there is a natural inner scale-i.e. the dissipation length-that can be 
used to non dimensionalize the equations. In this case, the Reynolds number is fixed, and 
of order unity, and the only remaining free parameter is the domain size. With this choice, 
increasing the resolution and increasing the domain size are entirely equivalent operations. 
It is then possible to define the characteristic scale of the solutions, once they are computed, 
in terms of, say, the (inverse) scale at which the velocity or magnetic spectrum peaks; the 
latter quantity being conceptually useful since it is related to the effective angular momen- 
tum transport. The problem of convergence then relates to the behavior of the characteristic 
scale of the solutions as the computational domain size increases. 

As we discussed in section 2, there are two limiting possibilities: one in which the 
solution scale diverges with the domain size, the other in which it becomes asymptotically 
independent of it. The former can lead to a net turbulent transport that is independent of 
diffusivity, the latter to a quasi-collisional transport that scales as some fixed multiple of the 
diffusivity. The evidence from the numerical simulations, in particular figure [U is that it is 
the latter and not the former that gives the correct asymptotic behavior. In other words, the 
characteristic scale of the solutions becomes independent of the size of the computational 
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Fig. 7. — Plot of the stresses as a function of Aa, the triangles represents the results of 
numerical calculations carried out with the MP5 scheme, with 128 points in the vertical 
direction. The horizontal dashed line represents the averaged stress value for the zero flux 
case; the dashed line shows a quadratic behavior; the solid line shows a linear behavior and 
the vertical solid line shows the linear stability boundary. 
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domain. If we accept this result as correct, there are a number of issues that naturally arise. 

The first concerns the physical interpretation of the result. In the weak field regime, 
MRI driven turbulence should be considered within the framework of a subcritical dynamo 
instability since the magnetic field necessary for the development of the MRI must be gen- 
erated self-consistently by the MRI turbulence itself. In dynamo theory, it is customary to 
distinguish between two types of dynamo action, small-scale and large-scale. Small-scale 
dynamos generate magnetic field with characteristic scale comparable to, or smaller than 
the characteristic scale of the velocity field. The only requirements for the operation of such 
dynamos are that the magnetic Reynolds number be sufficiently high, and that the under- 
lying velocity not be too symmetric. Clearly, both conditions can be satisfied here provided 
the system size is large. Large-scale dynamos, on the other hand, generate magnetic field 
with characteristic scale larger than that of the velocity. The latter types of dynamo are 
associated with flows lacking r eflectional symmetry , i.e. helical flows, and an inverse cascade 



of magnetic helicity. Recently, iTobias et al.l ( 120111 ) have argued that in the nonlinear regime 



it is often difficult to distinguish unambiguously between large- and small-scale dynamo ac- 
tion, and have introduced instead the idea of a system-scale dynamo, as one that produces 
magnetic structures comparable in size to the system scale irrespective of the scale of the 
velocity. The fact that, here, the scale at which the magnetic spectrum peaks is comparable 
with the dissipation scale and it is asymptotically independ ent of the system size indicates 



that the dynamo operating here is of the small-scale type (IVainshtein fc Kichatinov! 11986 



Boldyrev fc Cattaneoll2004j ). Dynamos of this type have no manifest inverse cascade of (mag- 
netic) helicity, they can lead to the production of magnetic energy but not of substantial 
magnetic flux. In the specific context of MRI driven turbulence a small-scale dynamo will 
not give rise to a "turbulent" transport of angular momentum. When regarded in the general 
scheme of dynamo theory, this result is actually not that surprising since the underlying flow 
is not strongly helical, and therefore, there is no a priori reason to expect an efficient inverse 
cascade. 

This brings us to the second issue, which relates to the astrophysical significance of these 
types of MRI solutions. Clearly, it is to be expected that this type of dynamo excitation 
will be prevalent in most (electrically conducting) discs. Its role in the transport of angu- 
lar momentum, and therefore in regulating the accretion rate, however will be negligible. 
With things as they stand, the only dynamo solutions that might have an impact on the 
angular momentum transport are those with an efficient inverse cascade. Interestingly, this 
implies that by their very nature such solutions should not be treated in the context of a 
local approximation. We find ourselves in a paradoxical situation in which if a solution is 
astrophysically relevant it cannot properly be described by a local model, and if it can be 
described by a local model then it is astrophysically irrelevant. The third issue concerns the 
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role of numerics in influencing what we understand, or think we understand in computational 
astrophysics. It is interesting that the asymptotic behavior that eventually emerges accord- 
ing to the present study is in many ways the most natural: that consonant with a subcritical 
small-scale dynamo instability developing in a system where there are not many reasons to 
expect otherwise. Likewise, the scaling behavior for the cases with weak imposed fluxes is 
asymptotically the one that makes most sense given the symmetries of the equations. Yet, 
their numerical realizations were somewhat long in the making, and had not been reported 
in the literature to date. The reason is that these asymptotic regimes requires a substantial 
dynamic range to manifest themselves. In this particular case, this was achieved by a com- 
bination of high resolution and high accuracy numerics. If we had tried to obtain the correct 
limiting behaviors with the less accurate codes, it would have been prohibitive^. Finally, it 
behooves us, to compare the conclusions of the present study with those of related studies in 
which different s c aling behaviors have been observed. Of particular relevance are the work of 



Kapyla &: Korpil ( 120101 ) who also consider unstrat ified shearing boxes but with non periodic 



boundary conditions; the work of lFromana (120101) who ha s the same set up as ours but with 
explicit diffusivities; and the work of iDavis et al.l (120101 ) who include stratification. In all 
these cases the authors report that the problem of convergence as stated above, is absent. 
In the present context, these conclusions can be attributed to one or more of the following 
possibilities: the simulations are not in the asymptotic regime because of limited resolution 
and/or accuracy; the extra physics has engendered an inverse cascade-i.e. the type of dy- 
namo action has changed from small-scale to large-scale; the extra physics has changed the 
nature of the subcritical dynamo instability and has caused it to become of the system-size 
type. The last two possibilities are related but distinct. They both refer to circumstances 
that lead to the production of large scale magnetic fields, but in a large-scale dynamo the 
velocity correlation length remains small, it is only the magnetic structures that become 
large; in a system-scale dynamo both the velocity and magnetic structur es grow in size unt il 
they occupy the largest available scale-the system size-hence the name (ITobias et al.ll201ll ). 



We begin by discussing the work of Kapyla fc Korpil ( 120101 ) . They consider unstratified 
shearing box simulations with boundary conditions at the upper and lower boundaries corre- 
sponding to impenetrable, stress-free velocities and vertical magnetic fields. The distinctive 
feature of these conditions is that they allow a flux of magnetic helicity in and out of the 
box. This has been advocated by several authors as a f eature that might be beneficial to 
large-scale dynamo action (see e.g. IVishniac fc Choi 120011 ) . The calculations are carried out 



1 It is useful to note that the scheme with quadratic reconstruction used here has order of accuracy 
comparable with ZEUS and ATHENA; two of the codes commonly used by the practitioners of numerical 
studies of the MRI 
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with explicit cliff usivities and the authors report that the effective angular momentum trans- 
port is substantial and insensitive to variations in the values of both the magnetic Reynolds 
and Prandtl numbers. The morphology of the solutions is dramatically different from the 
periodic cases, with two large current boundary layers forming in the vicinity of the upper 
boundaries and the generation of a coherent nearly uniform azimuthal field. All these findings 
are consistent with the activation of a large-scale type dynamo. There remains to be seen if 
these solutions persist at higher r esolutions, as we ll as h ow they map to a global geometry. 



These issues notwithstanding, the lKapyla fc Korpil (120101 ) solutions are quite remarkable and 



could potentially h ave considerable astrophysical significance. Next, we discuss stratification 



( IDavis et al.ll2010f ). Its presence drives two new effects: it introduces a new spatial scale in 
the problem, namely the vertical scale height that breaks the symmetry of the homogeneous, 
periodic models, and it also introduces buoyancy forces associated with thermal, pressure, 
and magnetic pressure fluctuations. It is not unlikely that these new effects could modify 
the type of dynamo action and drive flows and field coherent on scales comparable with the 
scale height. If this is indeed the case, it would be interesting to see how the solutions behave 
as the simulation size is increased for fixed scale height. Such a study would be numerically 
very demanding but astrophysically very useful. 



Finally, we consider the recent results of IFromangi ( I2010I ). They report that when ex- 



plicit diffusivities are included, the angular momentum transport becomes independent of 
resolution-i.e. the convergence problem goes away. If this result has a physical basis, it 
implies that the presence or absence of an inverse cascade and the type of associated dy- 
namo action (large-scale as as opposed to small-scale) depends on the detailed form of the 
diffusivities. The dynamo is small-scale for all three types of numerical diffusivities consid- 
ered here and in other similar studies, and large-scale for "physical" diffusivities. Although 
this possibility is not inconceivable, it is somewhat peculiar. It implies that the type of 
dynamo action is determined not by the symmetries of the problem, or large-scale features, 
but by the micro-physics. If correct it is a remarkable result. Another possibility, is that 
the simulatio ns are no t in th e asymptotic regime. To explore this possibility we note that 



the study by IFromangi (120101 ) was carried out by incorporating explicit diffusivities into the 
ZEUS code. In order for the explicit diffusivities to be correctly represented they must give 
rise to boundary layers whose size is not smaller than those arising from truncation errors. 
Thus the dynamic range of a code with numerical diffusivities gives an upper bound on the 
dynamic range achievable with the same code, the same resolution and explicit diffusivities. 
Now, ZEUS has an order of accuracy comparable to the code with quadratic reconstruc- 
tion, thus we can ask the following question. What is the resolution that we need with the 
quadratic code to see the same asymptotic behavior we observe with the MP5 code? This 
can be estimated by comparing the curves in Figure [1] and taking the ratio in resolutions 
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for which the curves corresponding to the PPM and MP5 code have similar slopes. This 
exercise gives a factor close to 4. A more careful analysis comparing boundary layers in a 
channel flow problem yields a similar estimate. Thus our conclusion is that our quadratic 
code requires a resolution exceeding 1000 grid-points to reproduce the asymptotic behavior 
observed by the MP5 code with 256 grid-points. Since the highest resolution in the work by 



Fromangi (120101 ) was 512 grid-points it is possible that their results may not yet be in the 



asymptotic regime. 
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